#### Replication archive - TESTING FOR NEGATIVE SPILLOVERS: IS PROMOTING HUMAN RIGHTS REALLY PART OF THE 'PROBLEM'?
#### Kelley, Simmons, Strezhnev
#### 

### Libraries
require(foreign)
require(ggplot2)
require(MASS)
require(Hmisc)
require(DataCombine)
require(ivpack)
require(ivmodel)
require(AER)
require(readstata13)
require(texreg)

### Helper Functions
# Clustered standard errors
cluster_se <- function(model, cluster){
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- length(model$coefficients) + length(model$zeta)
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

###### General 
set.seed(02138)

###################################
## Hafner-Burton (2008)
###################################

### Coverage data 
hb.coverage <- read.dta13("rrhrd9 small.dta")
hb.coverage <- subset(hb.coverage, !is.na(ccode))
hb.coverage <- hb.coverage[order(hb.coverage$ccode, hb.coverage$year),]

### Load in HB Shaming data
hb.data.raw <- read.dta13("shaming data.dta")

### Drop south Vietnam
hb.data.raw <- subset(hb.data.raw, idcode != 817)

### Load in CIRI data
ciri.data <- read.csv("ciri.csv")
### Fix missing data
ciri.data[ciri.data==-999] <- NA
ciri.data[ciri.data==-66] <- NA
ciri.data[ciri.data==-77] <- NA

### Subset to non-missing countries
hb.data <- subset(hb.data.raw, !is.na(idcode))
## Subset to only years before 2001 (shaming continues to 2000, but we have 2001's outcome data, so bonus obs)
hb.data <- subset(hb.data, year <= 2001&year >=1975)
### Order by country-year 
hb.data <- hb.data[order(hb.data$idcode, hb.data$year),]

### Reverse PIRI to be consistent
hb.data$PIRI <- 8 - hb.data$PIRI

########################
#### Country-specific fixes
########################

# Drop East Germany
hb.data <- subset(hb.data, idcode != 265)
# Fix 1991-2001 polity for Germany
hb.data[hb.data$idcode == 260&hb.data$year>=1991,]$Polity2 <- 10

#### Czech Republic 
hb.data[hb.data$idcode == 315&hb.data$year>=1993&hb.data$year<=2001,]$PIRI <- hb.data[hb.data$idcode == 316&hb.data$year>=1993&hb.data$year<=2001,]$PIRI 
hb.data <- subset(hb.data, idcode != 316)

### South Korea - missing coverage data
hb.data[hb.data$idcode == 732&hb.data$year<=2000&hb.data$year>=1975,]$Media_ai_rn <- hb.coverage[which(hb.coverage$ccode == 732),]$ainr
hb.data[hb.data$idcode == 732&hb.data$year<=2000&hb.data$year>=1975,]$Media_ai_br <- hb.coverage[which(hb.coverage$ccode == 732),]$aibr
hb.data[hb.data$idcode == 732&hb.data$year<=2000&hb.data$year>=1975,]$Media_ai_tot <- hb.coverage[which(hb.coverage$ccode == 732),]$aittl

### Lag year
hb.data <- slide(hb.data, Var="year", GroupVar="idcode", NewVar = "yearLag1", slideBy=-1)
### Verify we're not missing a year
hb.data[which(hb.data$year - hb.data$yearLag1 != 1),]

### Generate leads and lags of PIRI - lead 1 for equivalent of lagging remainder of obs
hb.data <- slide(hb.data, Var="PIRI", GroupVar="idcode", NewVar = "PIRILag1", slideBy=-1)
hb.data <- slide(hb.data, Var="PIRI", GroupVar="idcode", NewVar = "PIRILag2", slideBy=-2)
hb.data <- slide(hb.data, Var="PIRI", GroupVar="idcode", NewVar = "PIRILead1", slideBy=1)

hb.data <- slide(hb.data, Var="polrt", GroupVar="idcode", NewVar = "polrtLag1", slideBy=-1)
hb.data <- slide(hb.data, Var="polrt", GroupVar="idcode", NewVar = "polrtLag2", slideBy=-2)
hb.data <- slide(hb.data, Var="polrt", GroupVar="idcode", NewVar = "polrtLead1", slideBy=1)

### Keep 1981-2000 as the main data range
hb.data <- subset(hb.data, year >= 1981&year<=2000)

### Construct CIRI political rights index - association, speech

### Drop any without a CIRI score
hb.data.complete <- subset(hb.data, !is.na(PIRI))
### Safely drop the ones without a polrt score as well
hb.data.complete <- subset(hb.data.complete, !is.na(polrt))
#hb.data.complete <- subset(hb.data.complete, !is.na(POLRIGHT))

## Drop the ones without leads and lags
hb.data.complete <- subset(hb.data.complete, !is.na(PIRILag1))
hb.data.complete <- subset(hb.data.complete, !is.na(PIRILag2))
hb.data.complete <- subset(hb.data.complete, !is.na(PIRILead1))

hb.data.complete <- subset(hb.data.complete, !is.na(polrtLag1))
hb.data.complete <- subset(hb.data.complete, !is.na(polrtLag2))
hb.data.complete <- subset(hb.data.complete, !is.na(polrtLead1))

### Drop missing CAT, and ICCPR ratification
hb.data.complete <- subset(hb.data.complete, !is.na(cat_rat))
hb.data.complete <- subset(hb.data.complete, !is.na(ccpr_rat))

### Missing Polity2
hb.data.complete <- subset(hb.data.complete, !is.na(Polity2))

### Missing GDP PC
hb.data.complete <- subset(hb.data.complete, !is.na(gdppc))
hb.data.complete <- subset(hb.data.complete, !is.na(pop))
hb.data.complete <- subset(hb.data.complete, !is.na(Civilwar2))
hb.data.complete <- subset(hb.data.complete, !is.na(War))

### Binary democratic
hb.data.complete$democratic <- ifelse(hb.data.complete$Polity2 >= 6, 1, 0)

### Media criticism
hb.data.complete <- subset(hb.data.complete, !is.na(Media_ai_rn))

### Log version of the outcome?
hb.data.complete$Media_ai_tot_ph <- hb.data.complete$Media_ai_tot
hb.data.complete$Media_ai_tot_ph[hb.data.complete$Media_ai_tot_ph == 0] <- .01
hb.data.complete$logAI <- log(hb.data.complete$Media_ai_tot_ph) - mean(log(hb.data.complete$Media_ai_tot_ph))

### Did an NGO shame?
hb.data.complete$AIShame <- ifelse(hb.data.complete$Media_ai_rn > 0, 1, 0)

### Correlation between logged and binary measure
cor(hb.data.complete$logAI, hb.data.complete$AIShame)

### Sample size and number of countries
nCountry = length(unique(hb.data.complete$idcode))
sampleSize = nrow(hb.data.complete)

nCountry
sampleSize

###########################
#### Analysis
###########################

### Run the regression on PIRI - 3 lags
piri_reg_3 <- lm(PIRILead1 ~  AIShame + PIRI + PIRILag1 + PIRILag2 + cat_rat + ccpr_rat + democratic + log(gdppc) + log(pop) + Civilwar2 +  War + as.factor(year), data=hb.data.complete)
cluster_piri_reg_3 <- cluster_se(piri_reg_3, as.character(hb.data.complete$idcode))
coeftest(piri_reg_3, cluster_piri_reg_3)

### Run the regression on PIRI - 2 lags
piri_reg_2 <- lm(PIRILead1 ~  AIShame + PIRI + PIRILag1 + cat_rat + ccpr_rat + democratic + log(gdppc) + log(pop) + Civilwar2 +  War + as.factor(year), data=hb.data.complete)
cluster_piri_reg_2 <- cluster_se(piri_reg_2, hb.data.complete$idcode)
coeftest(piri_reg_2, cluster_piri_reg_2)

### Run the regression on PIRI - 1 lag
piri_reg_1 <- lm(PIRILead1 ~  AIShame + PIRI  + cat_rat + ccpr_rat + democratic + log(gdppc) + log(pop) + Civilwar2 +  War + as.factor(year), data=hb.data.complete)
cluster_piri_reg_1 <- cluster_se(piri_reg_1, hb.data.complete$idcode)
coeftest(piri_reg_1, cluster_piri_reg_1)

### Run the regression on polrt
polrt_reg <- lm(polrtLead1 ~  AIShame + polrt + cat_rat + ccpr_rat + democratic + log(gdppc) + log(pop) + Civilwar2 +  War + as.factor(year), data=hb.data.complete)
cluster_polrt_reg <- cluster_se(polrt_reg, hb.data.complete$idcode)
coeftest(polrt_reg, cluster_polrt_reg)

### Appendix Table 1
texreg(list(coeftest(piri_reg_3, cluster_piri_reg_3),coeftest(piri_reg_2, cluster_piri_reg_2),
            coeftest(piri_reg_1, cluster_piri_reg_1),coeftest(polrt_reg, cluster_polrt_reg) ), stars = c(0.01,  0.05, 0.10))

##### Plot mainline effects
##### Paper Figure 2
piri_effects <- data.frame(name = c("Physical Integrity\n3 Lags", "Physical Integrity\n2 Lags", "Physical Integrity\n1 Lag", "Political Rights\n1 Lag"), 
                           point = c(coeftest(piri_reg_3, cluster_piri_reg_3)["AIShame",1], coeftest(piri_reg_2, cluster_piri_reg_2)["AIShame",1], 
                                     coeftest(piri_reg_1, cluster_piri_reg_1)["AIShame",1], coeftest(polrt_reg, cluster_polrt_reg)["AIShame",1]),
                           se = c(coeftest(piri_reg_3, cluster_piri_reg_3)["AIShame",2], coeftest(piri_reg_2, cluster_piri_reg_2)["AIShame",2], 
                                  coeftest(piri_reg_1, cluster_piri_reg_1)["AIShame",2], coeftest(polrt_reg, cluster_polrt_reg)["AIShame",2])
)

piri_effects$lower95 = piri_effects$point - qnorm(.975)*piri_effects$se
piri_effects$upper95 = piri_effects$point + qnorm(.975)*piri_effects$se

piri_effects$name <- factor(piri_effects$name, levels=c("Physical Integrity\n3 Lags", "Physical Integrity\n2 Lags", "Physical Integrity\n1 Lag", "Political Rights\n1 Lag"))

p = ggplot(piri_effects, aes(y=point, x=name))
p = p + coord_flip()  
p = p + geom_hline(yintercept = 0,size=1,linetype="dashed") 
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 1)
p = p + scale_y_continuous(name="Estimated effect of NGO shaming on rights measure", limits=c(-.5, .5))
p = p + scale_x_discrete(name="")

theme_bw1 <- function(base_size = 20, base_family = "") {
  
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size*.9, colour = "black"),
      axis.text.y =       element_text(size = base_size, colour = "black"), # changes position of X axis text
      
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      plot.title =             element_text(face = "bold"),
      legend.position = "none"
    )
}
p = p + theme_bw()
print(p)
ggsave("figure2_hb_main_effects_replication.svg", width=8, height=3)
ggsave("figure2_hb_main_effects_replication.pdf", width=8, height=3)


#################################3
########### Piri with original logged shaming frequency measure
################################

### Run the regression on PIRI - 3 lags
piri_reg_3 <- lm(PIRILead1 ~  logAI + PIRI + PIRILag1 + PIRILag2 + cat_rat + ccpr_rat + democratic + log(gdppc) + log(pop) + Civilwar2 +  War + as.factor(year), data=hb.data.complete)
cluster_piri_reg_3 <- cluster_se(piri_reg_3, as.character(hb.data.complete$idcode))
coeftest(piri_reg_3, cluster_piri_reg_3)

### Run the regression on PIRI - 2 lags
piri_reg_2 <- lm(PIRILead1 ~  logAI + PIRI + PIRILag1 + cat_rat + ccpr_rat + democratic + log(gdppc) + log(pop) + Civilwar2 +  War + as.factor(year), data=hb.data.complete)
cluster_piri_reg_2 <- cluster_se(piri_reg_2, hb.data.complete$idcode)
coeftest(piri_reg_2, cluster_piri_reg_2)

### Run the regression on PIRI - 1 lag
piri_reg_1 <- lm(PIRILead1 ~  logAI + PIRI  + cat_rat + ccpr_rat + democratic + log(gdppc) + log(pop) + Civilwar2 +  War + as.factor(year), data=hb.data.complete)
cluster_piri_reg_1 <- cluster_se(piri_reg_1, hb.data.complete$idcode)
coeftest(piri_reg_1, cluster_piri_reg_1)

### Run the regression on polrt
polrt_reg <- lm(polrtLead1 ~  logAI + polrt + cat_rat + ccpr_rat + democratic + log(gdppc) + log(pop) + Civilwar2 +  War + as.factor(year), data=hb.data.complete)
cluster_polrt_reg <- cluster_se(polrt_reg, hb.data.complete$idcode)
coeftest(polrt_reg, cluster_polrt_reg)

## Results
texreg(list(coeftest(piri_reg_3, cluster_piri_reg_3),coeftest(piri_reg_2, cluster_piri_reg_2),
            coeftest(piri_reg_1, cluster_piri_reg_1),coeftest(polrt_reg, cluster_polrt_reg) ), stars = c(0.01,  0.05, 0.10))


##### Plot mainline effects
piri_effects <- data.frame(name = c("Physical Integrity\n3 Lags", "Physical Integrity\n2 Lags", "Physical Integrity\n1 Lag", "Political Rights\n1 Lag"), 
                           point = c(coeftest(piri_reg_3, cluster_piri_reg_3)["logAI",1], coeftest(piri_reg_2, cluster_piri_reg_2)["logAI",1], 
                                     coeftest(piri_reg_1, cluster_piri_reg_1)["logAI",1], coeftest(polrt_reg, cluster_polrt_reg)["logAI",1]),
                           se = c(coeftest(piri_reg_3, cluster_piri_reg_3)["logAI",2], coeftest(piri_reg_2, cluster_piri_reg_2)["logAI",2], 
                                  coeftest(piri_reg_1, cluster_piri_reg_1)["logAI",2], coeftest(polrt_reg, cluster_polrt_reg)["logAI",2])
)

piri_effects$lower95 = piri_effects$point - qnorm(.975)*piri_effects$se
piri_effects$upper95 = piri_effects$point + qnorm(.975)*piri_effects$se

piri_effects$name <- factor(piri_effects$name, levels=c("Physical Integrity\n3 Lags", "Physical Integrity\n2 Lags", "Physical Integrity\n1 Lag", "Political Rights\n1 Lag"))

p = ggplot(piri_effects, aes(y=point, x=name))
p = p + coord_flip()  
p = p + geom_hline(yintercept = 0,size=1,linetype="dashed") 
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_y_continuous(name="Estimated effect of 1 SD increase in NGO shaming on rights measure", limits=c(-.2, .2))
p = p + scale_x_discrete(name="")

theme_bw1 <- function(base_size = 20, base_family = "") {
  
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size*.9, colour = "black"),
      axis.text.y =       element_text(size = base_size, colour = "black"), # changes position of X axis text
      
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      plot.title =             element_text(face = "bold"),
      legend.position = "none"
    )
}
p = p + theme_bw()
print(p)
ggsave("extra_hb_main_effects_replication_log_scale.svg", width=8, height=3)
ggsave("extra_hb_main_effects_replication_log_scale.pdf", width=8, height=3)


###############
###  Demonstrate that the 2SLS results don't quite work 

hb.2sls <- ivreg(formula = PIRILead1 ~  polrtLead1 + PIRI + polrt  + cat_rat + ccpr_rat + democratic + 
                   log(gdppc) + log(pop) + Civilwar2 +  War + as.factor(year)|
                   logAI + PIRI  + polrt + cat_rat + ccpr_rat + democratic + 
                   log(gdppc) + log(pop) + Civilwar2 +  War + as.factor(year), data=hb.data.complete)

hb.2sls.cluster <- cluster.robust.se(hb.2sls, hb.data.complete$idcode)
hb.2sls.cluster

####################################
## Lupu replication
####################################

### Load in ICCPR data
lupu.data <- read.dta("ICCPR_FULL_DATA.dta")

### Load in CIRI data
ciri.data <- read.csv("ciri.csv")
### Fix missing data
ciri.data[ciri.data==-999] <- NA
ciri.data[ciri.data==-66] <- NA
ciri.data[ciri.data==-77] <- NA

### Merge with ICCPR data
lupu.data <- merge(lupu.data, ciri.data, by.x=c("ccode","year"), by.y=c("COW", "YEAR"), all.x=T, all.y=F)

### Sort by country code and year
lupu.data <- lupu.data[order(lupu.data$ccode, lupu.data$year),]

### Make outcomes factor variables
lupu.data$disapFactor <- factor(lupu.data$DISAP, ordered=T)
lupu.data$assnFactor <- factor(lupu.data$ASSN, ordered=T)
lupu.data$speechFactor <- factor(lupu.data$SPEECH, ordered=T)
lupu.data$new_relfreFactor <- factor(lupu.data$NEW_RELFRE, ordered=T)

### Make combined physical integrity/political rights


### 
### Lag each variable
lupu.data <- slide(lupu.data, Var="DISAP", GroupVar="ccode", slideBy=-1)
lupu.data <- slide(lupu.data, Var="ASSN", GroupVar="ccode", slideBy=-1)
lupu.data <- slide(lupu.data, Var="SPEECH", GroupVar="ccode", slideBy=-1)
lupu.data <- slide(lupu.data, Var="NEW_RELFRE", GroupVar="ccode", slideBy=-1)
lupu.data <- slide(lupu.data, Var="disapFactor", GroupVar="ccode", slideBy=-1)
lupu.data <- slide(lupu.data, Var="assnFactor", GroupVar="ccode", slideBy=-1)
lupu.data <- slide(lupu.data, Var="speechFactor", GroupVar="ccode", slideBy=-1)
lupu.data <- slide(lupu.data, Var="new_relfreFactor", GroupVar="ccode", slideBy=-1)

###
### Missing data - remove
lupu.data.final <- subset(lupu.data, !is.na(DISAP)) # Disappearances
lupu.data.final <- subset(lupu.data.final, !is.na(ASSN)) # Association
lupu.data.final <- subset(lupu.data.final, !is.na(SPEECH)) # Speech
lupu.data.final <- subset(lupu.data.final, !is.na(NEW_RELFRE)) # Religious Freedom

lupu.data.final <- subset(lupu.data.final, !is.na(iccpr_ratif)) # ICCPR
lupu.data.final <- subset(lupu.data.final, !is.na(injud)) # Injud
lupu.data.final <- subset(lupu.data.final, !is.na(polity)) # Polity
lupu.data.final <- subset(lupu.data.final, !is.na(durable)) # Durable

lupu.data.final <- subset(lupu.data.final, !is.na(civilwar)) # Civil War
lupu.data.final <- subset(lupu.data.final, !is.na(intlwar)) # International War
lupu.data.final <- subset(lupu.data.final, !is.na(gdpln)) # Log GDP
lupu.data.final <- subset(lupu.data.final, !is.na(popln)) # Log Pop
lupu.data.final <- subset(lupu.data.final, !is.na(INGO_uia)) # INGO UIA measure

lupu.data.final <- subset(lupu.data.final, !is.na(`DISAP-1`)) # Disappearances Lagged
lupu.data.final <- subset(lupu.data.final, !is.na(`ASSN-1`)) # Association Lagged
lupu.data.final <- subset(lupu.data.final, !is.na(`SPEECH-1`)) # Speech Lagged
lupu.data.final <- subset(lupu.data.final, !is.na(`NEW_RELFRE-1`)) # Religious Freedom Lagged

### Sample size
nrow(lupu.data.final) # Number of obs
length(unique(lupu.data.final$ccode)) # Number of countries

### Polr effects - Ordered logit

# Effect on Disappearances
ord_prob_Disap <- polr(disapFactor ~ `DISAP-1` + iccpr_ratif + injud + polity + durable + civilwar + intlwar + gdpln + 
                         popln + INGO_uia + as.factor(year), data=lupu.data.final, method="probit")
ord_prob_cluster_Disap <- cluster_se(ord_prob_Disap, lupu.data.final$ccode)

# Effect on Association
ord_prob_assn <- polr(assnFactor ~ `ASSN-1`  + iccpr_ratif + injud + polity + durable + civilwar + intlwar + gdpln + 
                        popln + INGO_uia + as.factor(year), data=lupu.data.final, method="probit")
ord_prob_cluster_assn <- cluster_se(ord_prob_assn, lupu.data.final$ccode)

# Effect on speech
ord_prob_speech <- polr(speechFactor  ~ `SPEECH-1`  + iccpr_ratif + injud + polity + durable + civilwar + intlwar + gdpln + 
                          popln + INGO_uia + as.factor(year), data=lupu.data.final, method="probit")
ord_prob_cluster_speech <- cluster_se(ord_prob_speech, lupu.data.final$ccode)

# Effect on Religious Freedom
ord_prob_relfre <- polr(new_relfreFactor  ~ `NEW_RELFRE-1`  + iccpr_ratif + injud + polity + durable + civilwar + intlwar + gdpln + 
                          popln + INGO_uia + as.factor(year), data=lupu.data.final, method="probit")
ord_prob_cluster_relfre <- cluster_se(ord_prob_relfre, lupu.data.final$ccode)

### Output to tex
texreg(list(coeftest(ord_prob_Disap, ord_prob_cluster_Disap),coeftest(ord_prob_assn, ord_prob_cluster_assn),
            coeftest(ord_prob_speech, ord_prob_cluster_speech),coeftest(ord_prob_relfre, ord_prob_cluster_relfre) ), stars = c(0.01,  0.05, 0.10))

### Linear model
### Disappearances
linear_disap <- lm(DISAP ~ `DISAP-1` + iccpr_ratif + injud + polity + durable + civilwar + intlwar + gdpln + 
                     popln + INGO_uia + as.factor(year), data=lupu.data.final)
cluster_linear_disap <- cluster_se(linear_disap , lupu.data.final$ccode)
coeftest(linear_disap, cluster_linear_disap)

# Effect on association
linear_assn <- lm(ASSN ~ `ASSN-1`  + iccpr_ratif + injud + polity + durable + civilwar + intlwar + gdpln + 
                    popln + INGO_uia + as.factor(year), data=lupu.data.final)
cluster_linear_assn <- cluster_se(linear_assn, lupu.data.final$ccode)
coeftest(linear_assn, cluster_linear_assn)

# Effect on speech
linear_speech <- lm(SPEECH ~ `SPEECH-1` + iccpr_ratif + injud + polity + durable + civilwar + intlwar + gdpln + 
                      popln + INGO_uia + as.factor(year), data=lupu.data.final)
cluster_linear_speech <- cluster_se(linear_speech, lupu.data.final$ccode)
coeftest(linear_speech, cluster_linear_speech)

# Effect on Religious Freedom
linear_relfre <- lm(NEW_RELFRE ~ `NEW_RELFRE-1` + iccpr_ratif + injud + polity + durable + civilwar + intlwar + gdpln + 
                      popln + INGO_uia + as.factor(year), data=lupu.data.final)
cluster_linear_relfre <- cluster_se(linear_relfre, lupu.data.final$ccode)
coeftest(linear_relfre, cluster_linear_relfre)

texreg(list(coeftest(linear_disap, cluster_linear_disap),coeftest(linear_assn, cluster_linear_assn),
            coeftest(linear_speech, cluster_linear_speech),coeftest(linear_relfre, cluster_linear_relfre)), stars = c(0.01,  0.05, 0.10))

#### Figure 3 - Lupu Replication
##### Plot mainline effects
main_effects <- data.frame(name = c("Disappearances", "Association", "Speech", "Religious Freedom"), 
                           point = c(coef(linear_disap)["iccpr_ratif"], coef(linear_assn)["iccpr_ratif"], coef(linear_speech)["iccpr_ratif"], coef(linear_relfre)["iccpr_ratif"]),
                           se = c(coeftest(linear_disap, cluster_linear_disap)["iccpr_ratif",2], coeftest(linear_assn, cluster_linear_assn)["iccpr_ratif",2],
                                  coeftest(linear_speech, cluster_linear_speech)["iccpr_ratif",2], coeftest(linear_relfre, cluster_linear_relfre)["iccpr_ratif",2])
)

main_effects$lower95 = main_effects$point - qnorm(.975)*main_effects$se
main_effects$upper95 = main_effects$point + qnorm(.975)*main_effects$se

p = ggplot(main_effects, aes(y=point, x=name))
p = p + coord_flip()  
p = p + geom_hline(yintercept = 0,size=1,linetype="dashed") 
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 1)
p = p + scale_y_continuous(name="Estimated effect of ICCPR ratification on average CIRI rights score", limits=c(-.2, .2))
p = p + scale_x_discrete(name="")

theme_bw1 <- function(base_size = 20, base_family = "") {
  
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size*.9, colour = "black"),
      axis.text.y =       element_text(size = base_size, colour = "black"), # changes position of X axis text
      
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      plot.title =             element_text(face = "bold"),
      legend.position = "none"
    )
}
p = p + theme_bw()
print(p)
ggsave("figure3_lupu_main_effects_replication.svg", width=8, height=3)
ggsave("figure3_lupu_main_effects_replication.pdf", width=8, height=3)

###############################
### Is there a trade-off between association and disappearances - sensitivity Analysis
####################################

### First stage regression under assumption of no direct effect
lupu.2sls <- ivreg(formula = DISAP ~ ASSN + injud + polity + durable + civilwar + intlwar + gdpln + 
                     popln + INGO_uia + `DISAP-1` + `ASSN-1` + as.factor(year) |
                     iccpr_ratif + injud + polity + durable + civilwar + intlwar + gdpln + 
                     popln + INGO_uia + `DISAP-1` + `ASSN-1` + as.factor(year),
                   data = lupu.data.final)
# Clustered SEs
lupu.2sls.cluster <- cluster.robust.se(lupu.2sls, lupu.data.final$ccode)

# Lupu first-stage regression
lupu.first.stage <- lm(ASSN ~ iccpr_ratif + injud + polity + durable + civilwar + intlwar + gdpln + 
                         popln + INGO_uia + `DISAP-1` + `ASSN-1` + as.factor(year), data=lupu.data.final)
cluster_fs <- cluster_se(lupu.first.stage, lupu.data.final$ccode)
coeftest(lupu.first.stage,cluster_fs)

# Get a sense of the magnitude of the first-stage (standardize by pooled sd)
sd_pooled_assn <- sqrt(((sum(lupu.data.final$iccpr_ratif == 1) - 1)*var(lupu.data.final$ASSN[lupu.data.final$iccpr_ratif == 1]) + 
                          (sum(lupu.data.final$iccpr_ratif == 0) - 1)*var(lupu.data.final$ASSN[lupu.data.final$iccpr_ratif == 0]))/(nrow(lupu.data.final) - 2))

stdEffect.fs <- coeftest(lupu.first.stage,cluster_fs )["iccpr_ratif",1]/sd_pooled_assn

### Sensitivity parameters Cohen's d
cohens_d <- c(seq(-.25, .25, by =.01), stdEffect.fs)
## Standard deviation
sd_pooled <- sqrt(((sum(lupu.data.final$iccpr_ratif == 1) - 1)*var(lupu.data.final$DISAP[lupu.data.final$iccpr_ratif == 1]) + 
                     (sum(lupu.data.final$iccpr_ratif == 0) - 1)*var(lupu.data.final$DISAP[lupu.data.final$iccpr_ratif == 0]))/(nrow(lupu.data.final) - 2))

### What's the magnitude of the effect on ASSN
coeftest(linear_assn, cluster_linear_assn)["iccpr_ratif",1]

#################
## Sensitivity analysis figure
#################

### Store estimates
point_result <- rep(NA, length(cohens_d))
se_result <- rep(NA, length(cohens_d))

### What data
working.data <- lupu.data.final

## for each sensitivity parameter value
for (sens in 1:length(cohens_d)){
  
  ### Effect magnitude
  DE <- cohens_d[sens]*sd_pooled
  
  ### Subtract off the direct treatment effect
  working.data$DISAPDeBiased <- working.data$DISAP - DE*working.data$iccpr_ratif
  
  ### Instrumental regression on de-biased data
  boot.2sls <- ivreg(formula = DISAPDeBiased ~ ASSN + injud + polity + durable + civilwar + intlwar + gdpln + 
                       popln + INGO_uia + `DISAP-1` + `ASSN-1` +  as.factor(year) |
                       iccpr_ratif + injud + polity + durable + civilwar + intlwar + gdpln + 
                       popln + INGO_uia + `DISAP-1` + `ASSN-1` + as.factor(year),
                     data = working.data)
  
  
  boot.cluster <- cluster.robust.se(boot.2sls, working.data$ccode)
  ### Save effect
  point_result[sens] <- boot.cluster["ASSN", 1]
  se_result[sens] <- boot.cluster["ASSN", 2]
  
}

### Sensitivity analysis plot - Figure 4
plot.frame <- data.frame(point=point_result, se=se_result, sensitivity=cohens_d)
plot.frame$lower95 <- plot.frame$point - qnorm(.975)*plot.frame$se
plot.frame$upper95 <- plot.frame$point + qnorm(.975)*plot.frame$se

p <- ggplot(data=plot.frame, aes(x=sensitivity, y=point)) + 
  geom_line( lwd= 1) +
  geom_line(aes(y=upper95), linetype=2, color="grey50") +
  geom_line(aes(y=lower95), linetype=2, color="grey50") +
  geom_hline(yintercept=0, linetype=3, lwd=1.5,  color="black") +
  geom_vline(xintercept=0,  linetype=3, lwd =1.5, color="black") +
  geom_segment(y = plot.frame$upper95[plot.frame$sensitivity == stdEffect.fs], x=stdEffect.fs, yend=plot.frame$lower95[plot.frame$sensitivity == stdEffect.fs], xend=stdEffect.fs, color="darkgrey", lwd=1.5, lty=1) + 
  geom_segment(y = plot.frame$upper95[plot.frame$sensitivity == 0], x=0, yend=plot.frame$lower95[plot.frame$sensitivity == 0], xend=0, color="lightgrey", lwd=1.5, lty=1) +
  
  scale_x_continuous("Standardized direct effect of ICCPR Ratification on CIRI Disappearance") +
  scale_y_continuous("Effect of CIRI Association on CIRI Disappearance") +
  
  theme_bw()
print(p)
ggsave("figure4_trade_off_plot_assn_disappearance_lupu.svg", width=8, height=5.5)
ggsave("figure4_trade_off_plot_assn_disappearance_lupu.pdf", width=8, height=5.5)


